#include <fstream>
#include <istream>
#include <sstream>
#include <vector>
#include <deque>
#include <string>
#include <iostream>



#define DEBUG_LEVEL 6
#include <debug_func.h>

#include <seq_pairs.h>
#include <blast_results.h>
#include <bcpp_read_seq.h>
#include <codeml_results.h>
#include <progress_indicator.h>
#define _MY_STL_MACROS
#include <stlhelper.h>
#include <ultoa.h>
#include <print_error.h>

#include <boost/foreach.hpp>

#include "get_arg.h"


using std::ifstream;
using std::istream;
using std::istringstream;
using std::string;
using std::cerr;
using std::cin;
using std::cout;
using std::runtime_error;
using BioCPP::t_collect_sequences;
using BioCPP::t_collect_sequences_report_errors_to_strm;
using BioCPP::eSEQ_NODUPLICATES;
using BioCPP::eSEQ_DUPLICATES;
using bcpp_align::t_align_results;


//________________________________________________________________________________________

//		do_parse_kaks

//________________________________________________________________________________________
struct do_parse_kaks
{
	t_collect_sequences&					seq_db;
	t_program_args&							args;
	std::string								comment_letters;
	do_parse_kaks(	t_collect_sequences&	sequences_,
					t_program_args&			args_,
					const std::string&		comment_letters_):
					seq_db(sequences_),
					args(args_),
					comment_letters(comment_letters_){}
	void operator()(istream& is, unsigned cnt_line)
	{
		// get data for one protein pair
		string name1, name2;

		// comment if false
		string comment_line;
		if (!seq_pairs::parse_ids(is, name1, name2, comment_line, comment_letters,
								  seq_pairs::ePARSE_IDS_ONLY))
			return;

		t_codeml_results dndsdata;
		is >> dndsdata;
		t_align_results aligndata;
		is >> aligndata;

		if (!is)
			throw runtime_error("The kaks distances data is not in the correct format");

		// get protein and gene indices
		unsigned seq_id1, seq_id2;
		if (!seq_db.get_index_of_name(name1, seq_id1))
			throw runtime_error("The sequence for ["  + name1 + "] is missing.");
		if (!seq_db.get_index_of_name(name2, seq_id2))
			throw runtime_error("The sequence for ["  + name2 + "] is missing.");

		string seq1 = seq_db.sequence(seq_id1);
		string seq2 = seq_db.sequence(seq_id2);
		insert_gaps_into_sequence_pair(aligndata, seq1, seq2,
										false, name1, name2);
		unsigned identical = 0;
		for (unsigned i = 0; i < seq1.size(); ++i)
			if (seq1[i] == seq2[i])
				++identical;
		aligndata.percent_id = identical * 1.0 / seq1.length();
		cout
			<< name1 << "\t" << name2 << "\t" 
			<< dndsdata << "\t"
			<< aligndata << "\n";
	}
};






//________________________________________________________________________________________

//		read_sequences

//________________________________________________________________________________________
void read_sequences(t_collect_sequences& seq_db, t_program_args& args)
{
	{
		ifstream f_sequences(args.file_seq.c_str());
		if (!f_sequences)
			throw runtime_error("Could not open " + args.file_seq +
										" to retrieve cDNA sequences");

		args.VERBOSE()	<< "\tReading protein sequences  from [ "
						<< args.file_seq << " ]...\n";
		t_progress_indicator dots(args.VERBOSE(), 1000);
		t_collect_sequences_report_errors_to_strm report_errors(args.CERR());
		seq_db.init_from_stream(f_sequences, ".*", dots, report_errors);
	}
	args.VERBOSE() << "\tDone...\n";
	args.VERBOSE() << "\t" << seq_db.size() << " sequences\n";
}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	t_program_args args;

	try
	{

		debug_func_cerr(5);

		// don't use C streams
		std::ios::sync_with_stdio(false);


		// get command line arguments
		if (!process_command_line_options(argc, argv, args))
			return 1;


		//
		// read sequences
		//
		t_collect_sequences seq_db;
		read_sequences(seq_db, args);

		do_parse_kaks parse_kaks(seq_db, args, args.comment_letters);
		seq_pairs::read_data(parse_kaks, args.file_input, args.VERBOSE());
	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}

